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Abstract 

A tutorial introduction to the technique of Molecular Dynamics (MD) 
is given, and some characteristic examples of applications are described. 
The purpose and scope of these simulations and the relation to other sim- 
ulation methods is discussed, and the basic MD algorithms are described. 
The sampling of intensive variables (temperature T, pressure p) in runs 
carried out in the microcanonical (NVE) ensemble (7V= particle number, 
V = volume, E — energy) is discussed, as well as the realization of other 
ensembles (e.g. the NVT ensemble). For a typical application example, 
molten Si02, the estimation of various transport coefficients (self-diffusion 
constants, viscosity, thermal conductivity) is discussed. As an example of 
Non-Equilibrium Molecular Dynamics (NEMD), a study of a glass-forming 
polymer melt under shear is mentioned. 
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1 Introduction: Scope and Purpose of Molecu- 
lar Dynamics Simulations 



1.1 Molecular Dynamics and its Relation to Other Meth- 
ods of Computer Simulation 

Computer simulations in condensed matter physics aim to calculate structure 
and dynamics from atomistic input . The theoretical basis of this ap- 

proach is statistical thermodynamics. The conceptually simplest approach is 
the classical Molecular Dynamics (MD) method 'W'W'T: one simply solves nu- 
merically Newton's equations of motion for the interacting many particle system 
(atoms or molecules interacting, e.g., with pair potentials). The basics of the 
method thus is nothing but classical mechanics, and one creates a deterministic 
trajectory in the phase space of the system. Now the idea is to simply take 
time averages of the observables of interest along this trajectory, and rely on 
the ergodicity hypothesis of statistical mechanics, which asserts that these time 
averages are equivalent to ensemble averages of the appropriate microcanonical 
(NVE) ensemble. Of course, Newton's equations of motion conserve the to- 
tal energy E, and hence the conjugate intensive thermodynamic variables such 
as temperature T and pressure p can only be inferred indirectly and exhibit 
fluctuations (since the particle number N is finite and sometimes even fairly 
small, such fluctuations must not be neglected and need careful consideration). 
Sometimes it is advantageous to directly realize other ensembles of statistical 
mechanics, such as the constant volume V- constant temperature T (NVT) 
ensemble or the NpT ensemble, and — as we shall see later — this is indeed 
feasible by introducing a coupling to appropriate "thermostats" or "barostats" . 

An alternative way of carrying out a MD simulation at constant tempera- 
ture is possible by introducing an artificial weak friction force, together with 
random forces whose strength are controlled by the fluctuation-dissipation the- 
orem. Such techniques are very common e.g. for the simulation of polymer 
melts IHIE]- This method is very closely related in spirit to stochastic simula- 
tion methods such as "Brownian Dynamics" where one simulates a Langevin 
equation (the inertial term in the equation of motion being omitted). While 
dynamical correlations for such methods differ somewhat from strictly micro- 
canonical MD methods, for the computation of static properties from averages 
along the stochastic trajectory in phase space such methods can be advanta- 
geous. 

This statement is also true for the importance sampling Monte Carlo (MC) 
method As is well known MC sampling means that one creates 

a random walk-like trajectory in configuration space, controlled by transition 
probabilities that ensure the approach to thermal equilibrium via a detailed 
balance condition. Many of the practical aspects of computer simulations, such 
as "statistical errors" and systematic errors due to the finite size of the simulated 
system or the finite "length" of the simulated trajectory (or observation time, 
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respectively), are shared by all these simulation methods. 

However, one important caveat needs to be made: it is quantum mechanics 
that describes the basic physics of condensed matter, and not classical mechan- 
ics. However, attempting a numerical solution of the Schrodinger equation for 
a system of many nuclei and electrons still is premature and not at all feasible 
even at the fastest computers. Thus, one has to resort to approximations. One 
very popular approach is the "ab initio MD" or "Car-Parrinello- method" |12| . 
where one includes some electronic degrees of freedom into MD via density 
functional theory (DFT) ^21 • The huge advantage of this technique is that one 
no longer relies on effective interatomic potentials, which often are only phe- 
nomenologically chosen ad hoc assumptions, lacking any firm quantum-chemical 
foundation. However, the disadvantage of this technique is that it is several or- 
ders of magnitude slower than classical MD, and hence only very short time 
scales and very small systems are accessible. Furthermore, the method is un- 
suitable to treat systems with van der Waals like forces, such as in rare gases, 
where one still is better off with the simple Lennard- Jones potential (perhaps 
amended by three-body forces) pQ. Also, normally ionic degrees of freedom 
are still treated classically. Alternatively, one can still use effective potentials 
between ions and/or neutral atoms as in classical MD or MC, but rely on quan- 
tum statistical mechanics for the ionic degrees of freedom: this is achieved by 
path integral Monte Carlo (PIMC) |14II15| or path integral molecular dynamics 
(PIMD) ^lEIEl- Such techniques are indeed crucial for a study of solids at 
low temperatures, to ensure that their thermal properties are compatible with 
the third law of thermodynamics. For most fluids (of course, quantum liquids 
such as ^He and ^He are an exception) classical MD is fine, and will henceforth 
be considered exclusively in this article. 

What information do we then desire to extract from the simulations? When 
we consider systems in thermal equilibrium, the first task is to calculate static 
properties. E.g., in a fluid a basic property is the static structure factor S{k) jJU] 

5(fc) = (|<5p(fc)p}, (1) 

where 5p{k) is a spatial Fourier transform of density fluctuations, k being the 
wavevector. In addition, one wants to calculate time-dependent correlation 
functions, that describe the decay of small thermal fluctuations with time. A 
quantity of this type that will be discussed later is the intermediate scattering 
function S{k, t), 

S{k,t) = {Sp{-k,0)Sp{k,t)). (2) 

If one considers systems out of thermal equilibrium, an important applica- 
tion of MD is the study of systems exhibiting a steady state flow under shear 
deformation. The purpose of such NEMD [20112 1| work can be the estimation 
of transport coefficients (e.g. the shear viscosity) if the deformation is weak 
enough so that the system is still in the regime of linear response jj^l- However, 
also the study of nonlinear phenomena (such as "shear thinning", a decrease 
of the effective viscosity with increasing shear rate [2011^) is of interest j22| . 
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In addition, one can also study non-steady state transient behavior, as it oc- 
curs e.g. after a sudden quench from one state to another, where one wishes 
to study the approach of the system to its new thermal equilibrium. Classical 
examples of this problem are nucleation of fluid droplets in a supersaturated 
gas or the kinetics of phase separation in a binary mixture ("spinodal decom- 
position" However, often such processes are too slow, and cannot really 
be studied by NEMD, and instead one has to use simulations of coarse-grained 
models by Non-Equilibrium Monte Carlo (NEMC) piUI^ITT]. 

Now this list of simulation techniques and problems that one may simu- 
late sounds wonderful - everything can be studied with computer simulation 
methods, even problems outside the standard realm of physics: the sponta- 
neous formation of traffic jams on highways |24j . anomalous time-dependent 
autocorrelation functions of heart beats of human beings suffering from heart 
diseases [H], critical fluctuations at the stock market before a crash or 
shock waves that can destroy a silo used in agriculture when the grains of corn 
stored in it flow out at the wrong speed |27], etc. All these problems are in fact 
simulated by theoretical physicists using techniques very similar to MD or MC. 
However, one must nevertheless always keep in mind that computer simulations 
- as all other techniques! - suffer from some very important technical limita- 
tions that the practitioner always must be aware of, just as an experimentalist 
in an inelastic scattering experiment must be aware of that the resolution of his 
energy and momentum transfer measurements limits his data analysis, and that 
in addition there are statistical errors due to limited intensity of the radiation 
source. 

1.2 Limitations of Computer Simulations: 

When to Choose Which Method? 

The main limitation of atomistic simulations comes from the fact that one often 
must bridge many decades in spatial scale and even more decades in time scale 
to describe the phenomena of interest 28..29| . As an example, we discuss here 
the problem of phase separation of a mixture of a polymer (such as polybu- 
tadiene) and its deuterated counterpart |30| . If one carries out a quenching 
experiment, suddenly reducing the temperature from a high value in the one 
phase region to a lower value inside the miscibility gap and records the equal- 
time structure factor S{k), Eq. (^), at different times t after the quench, one 
observes the growth of a peak at a position km (t) [30] . Now polybutadiene is a 
flexible linear macromolecule, which in the molten state forms random walk-like 
coils that exhibit nontrivial structure from the scale of covalent C-C and C-H 
bonds (i.e., of the order of 1 A) to the coil radius (which is of the order of 10^ A, 
for the molecular weights used in the experiment). The collective length scale 
i{t) « 2n /km.{t) is of the order of 1000 A already in the initial stage of phase 
separation, however. Clearly, such a study of cooperative phenomena of a large 
number of chains would be prohibitively difficult, if we would try a chemically 
realistic, atomistically detailed description. Moreover, the description of effec- 
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tive potentials driving this phase separation between protonated and deuterated 
polymers which otherwise are chemically identical is quite subtle: in the frame- 
work of classical statistical mechanics, the masses of the particles cancel out 
from all static averages, and hence such a mixture would be an ideal mixture, 
perfectly miscible, no phase separation would ever occur. The reason that phase 
separation is observed in the experiment [30] is a quantum effect, the zero point 
vibrational motions of hydrogens and deuterons differ: the lighter hydrogens 
need more space and this causes an effective repulsive interaction between un- 
like species. 

An even greater problem is the disparity of time scales encountered in this 
example: while the time-scale of bond length and bond angle vibrations is of the 
order of 10^^'^ s, already the time needed for a coil to renew its configuration can 
be estimated as 10~^ s, 8 orders of magnitude larger than the vibration times, 
for the considered molecular weight and temperature of the experiment. 
The collective dynamics, however, is even much slower, because the thermody- 
namic driving forces are very weak. Thus, the experiment shows |3(J| that the 
interesting phenomena happen on the time scales from a few seconds to 1000 
seconds, when a scattering peak develops at km{t) and first grows more or less 
at the same position and then shifts to smaller and smaller wavevectors as the 
coarsening of the structure proceeds. And while for the case of phases separa- 
tion in metallic alloys the situation is better with respect to the length scale 
km{t) |23I31| . km{t) ~ 0.01 — 0.1 A~^, the typical time scale is still from 0.1 s to 
10^ s, and hence for a chemically realistic molecular dynamics simulation, with 
a time step 6t in the range from 1 fs to 10 fs, the task is quite hopeless, one 
would have to simulate over a range of 10^^ time steps which is many orders of 
magnitude more than what is feasible nowadays. 

Such slow phenomena as spinodal decomposition in solid metallic alloys or 
fluid polymer mixtures can only be simulated by very simplified coarse-grained 
models, where chemical detail is sacrificed, and one applies non-equilibrium 
Monte Carlo (NEMC) methods rather than MD. These coarse-grained simula- 
tions nevertheless are very useful, both for solid alloys 32 and for polymers |33| . 
In the latter case, several successive chemical monomers along the chain are in- 
tegrated into one "effective bond" . Moreover, one also simulates relatively short 
(unentangled) chains, and hence one does not attempt to study chains of large 
molecular weight as done in the experiment |30j . Even with these simplifica- 
tions, it is difficult to deal with such long wavelength phenomena: the computer 
simulation can never deal directly with the system in the thermodynamic limit, 
one always can treat only a finite system. In this case of spinodal decompo- 
sition of binary polymer mixtures a cubic box containing a few hundred 
or a few thousand short polymer chains is studied. To avoid surface effects at 
the boundaries of the simulation box, periodic boundary conditions are used. 
Such systems sometimes are called quasi-infinite, but this is not quite true: the 
finite size of the system still has notable effects, e.g. the reciprocal space is 
"quantized" , since the only wavevectors that are compatible with the periodic 
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boundary conditions are of the form 



) = -j-{nx,ny,nz) 



(3) 



where L is the size of the box and nx,ny, Uz are integers. 

In practice one does find for this problem that the position km {t) where the 
peak grows occurs near values of k where /i < 10, and hence the discreteness of 
fc-space is a real practical problem. Despite such problems, the simulations of 
collective phenomena are useful, but it is always advisable to carefully 

pay attention to possible artifacts caused by the finite size of the simulation 
box. 

1.3 Internal Dynamics of Polymer Chains: An Example 
on What Can be Done When Models are Carefully 



This example of unmixing in polymer blends |33j is not meant to imply that 
MD simulations for polymers are not feasible at all: on the contrary, MD work 
for polymers can be very useful and can even be compared to experiment quan- 
titatively, but only for carefully selected problems! This is shown in an exam- 
ple |34[I35| dealing with the relaxation of the configuration of polymer coils. 
While for entangled chains with a degree of polymerization of the order of 10'' 
the relaxation time of the coil configuration typically is at least 10~^ s, the 
choice of shorter non-entangled chains such as Cioo H202 brings down to 
about Tfl « 1 ns, if a sufficiently high temperature is chosen (T = 509 K in 
this case). Fig. shows a comparison of data for the single-chain (normalized) 
coherent intermediate scattering function 



obtained by the neutron spin echo technique, with MD simulation results of 
a suitable model. In Eq. |@}, g is the scattering vector, Np is the degree of 
polymerization of the chain {Np — 100 here), and fm{to) is the position of the 
m'th monomer of a chain (1 < m < Np) at time to. The average (• ■ ■ ) in the 
simulation is taken over all chains in the system and is also meant to represent 
the time average over the time to. Note that we have already made use of the 
fact that in thermal equilibrium a time-displaced correlation function as written 
on the rhs of Eq. Q depends only on the difference t between the two times 
t + tg, to, and not on these two times separately. 

It is seen from Fig.nthat there is an almost perfect agreement between the 
experimental data and the simulations, over two decades in time and almost a 
decade in the wavevector, and in this scaled plot there is no adjustable param- 
eter whatsoever! Thus, this agreement is by no means a trivial fitting exercise. 



Chosen 
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but rather it is truly significant. While the original conclusion of the experimen- 
talists was that their data prove the validity of the Rouse model |36[l37j . one 
now knows, thanks to the simulation (Fig.|2Jl, that this is not quite correct. Re- 
member that the Rouse model describes the polymer as a harmonic chain, which 
experiences friction and stochastic forces which represent the interactions with 
the surrounding chains. This simplistic model contains only two parameters, an 
effective size cr of a bead, and the friction coefficient C- Both are known from 
independent estimates of other quantities in the simulation (static structure fac- 
tor S{q, 0) and self-diffusion constant D^j of the chains, respectively). So again 
one can perform a comparison between the simulation and the theory without 
any adjustable parameters whatsoever (Fig.[2Jl. One sees that the Rouse model 
works very well for small values of q, while it fails for larger wavelengths. Later 
experimental work did confirm the conclusion of the simulations p?5] that the 
Rouse model is indeed not so perfect as originally claimed. 

This figure also shows that S{q, t) for typical values of q really does decay 
to zero on the scale of 1 to 10 nanoseconds. We emphasize again that this 
was possible only due to carefully chosen simulation parameters: C100H202 is 
a relatively short chain, and T = 509 K is a relatively high temperature. The 
experiment deliberately was done for such parameters to allow a meaningful 
comparison with a simulation. As was said above, for other parameters it could 
easily happen that the relaxation times for q « i?^ ^ {Rg being the gyration 
radius of the chains) is larger by a factor of 10^-10 than in the present case. 
While at large enough q it still may be possible to study the relaxation in 
the "time window" accessible for the scattering experiment, it would not make 
sense to compare with a computer simulation: the latter cannot reach thermal 
equilibrium if the Rouse time (the time needed to equilibrate the configuration 
of the whole chain) is so large. Note that often it is thought that there is always a 
direct correspondence between quantities that one can study by inelastic neutron 
scattering and the corresponding MD observations, since both methods have a 
"time window" of about 1 ns. However, in an experiment one can invest a time 
of days or even weeks to carefully anneal the sample and thus prepare very good 
equilibrium. But this is not the case in a simulation: in MD work one almost 
always has to start from some initial state which is not yet characteristic for the 
desired equilibrium, and let the system relax towards equilibrium with a MD 
run that also can last only a few nanoseconds! Therefore meaningful inelastic 
scattering experiments can study the relaxation of fast degrees of freedom in 
fluids rather close to the glass transition temperature of polymers or other glass- 
forming systems, since the system is in good thermal equilibrium. In contrast 
to this MD simulations of a corresponding model can reach equilibrium only 
at much higher temperatures, however. These simple considerations are almost 
obvious but nevertheless often ignored - therefore we emphasize these points 
here. 

A related misconception is that a MD simulation is the better the more 
chemical detail is included in the model: however, this attitude is completely 
wrong! In fact, the level of detail in the simulations of Paul et al. |34II35| 
normally did neither include fluctuations in the lengths of the C-C bonds along 
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the chain backbone, nor were the hydrogen atoms exphcitly included: the bond 
length was kept at the experimental value ice = 1-53 A, and the CH2 groups 
were treated as "united atoms" (also no full distinction between the CH3 groups 
at the chain ends and the CH2 groups in the interior was made). If one includes 
the hydrogens explicitly, the program is about a factor 10 slower, and a modest 
gain in accuracy of the model is more than outweighed by about 3 times larger 
statistical errors. Note that this dramatic slowing down of the code is inevitable 
due to the very stiff potentials that need to be included in such an "all atom"- 
calculation and which necessitate then a particularly small time step (and also 
the number of atoms is 3 times larger). The potentials actually used for the 
simulations shown in Figs. and [21 are bond angle potentials U{Q) and torsion 
potentials f/(^!>) of the form 

1 ^ 
U{e) = e - cos eo)2 , t/(</)) = ^ a„ cos"(</)) , (5) 

n=0 

while non-bonded monomers interact with a Lennard-Jones potential 

U{n,)^Aec.fj[ia/nj)'^ ~ia/n,f] , a,/3e {CH2,CH3} . (6) 

The parameters fee, Qo, o-n, £ap, c are given in "3T'35 . In principle, such effec- 
tive potentials for classical MD simulations should be deduced from quantum- 
chemical calculations. However, polyethylene is a much too large molecule to 
do that: only the bond angle and torsional potentials have a quantum chemical 
basis, although even on this issue there is no full consensus between differ- 
ent groups. Of course, it is clear that the Lennard-Jones potential used here, 
{Eq. (|SJ)} is completely ad hoc, and the parameters are just optimized in order 
to fit as many experimental data as possible. For the case of C100H202, a- box 
of linear dimension i = 50 A allows to include 40 chains in the simulation, and 
finite size effects are then negligible for the wavevectors of interest. It is clear 
that a treatment of longer chains would not only require larger box sizes but 
also much longer runs, and therefore are very difficult. We also emphasize that 
Eq. (jSJ is not a good basis to describe interactions between chemically dissimilar 
polymers with sufficient accuracy. 



2 MD algorithms and some simulation "knowhow" 

2.1 Classical Mechanics and the Verlet Algorithm 

We consider a system of N particles (atoms) with Cartesian coordinates X = 
{fi}, « = 1, • • • , TV, in a d-dimensional space. The dynamics then is described 
by Newton's equations of motion, 

■m^Ti = —r- = Ji , (7) 

o n 
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mi being the mass of the i'th particle, and fi the force acting on it, which 
we assume to be entirely due to interactions with other particles. Thus, the 
potential Upot(X) = [/pot(n., • ■ ■ ,nv) is written as 

N-l N 

where in the last step we have furthermore made the simplifying assumption 
that C/pot is pairwise additive (this assumption, which is reasonable for simple 
liquids, is not really necessary, it can be generalized whenever appropriate, we 
only make it for the sake of simplicity of our presentation). Thus 

= - E dU{n,)/dr^, = ■ (9) 

The total energy 



^ 1 



i=l 

is a constant of motion, 



dE ^ • •• ^ . ^ 

— = ^ nun ri-^U- fi = ■ (11) 

1=1 i=l 



MD now means that Newton's equations of motion are integrated numeri- 
cally. A computationally efficient scheme is the so-called Verlet algorithm ■ 
To derive it, we expand •ri{t) forward and backward in time. 



n{t + St) = n{t) + Sm{t) + i^{5tfm + \{5t)%{t) + 0{{5tf), (12) 



n{t - St) ^ nit) - Stnit) + T^{St)^U{t) - ]:{St)%{t) + 0{{St)^) . (13) 

Here Vi{t) is the velocity of the i'th particle at time t, and hi{t) is the corre- 
sponding vector that appears in the third order of the Taylor expansion. 

For r{t) we have already substituted Newton's equation. Adding both ex- 
pansions, the odd terms cancel, and hence 

r,{t + St) = 2nit) - n{t - St) + —{St)^f:,{t) + OiiSt)^). (14) 

nil 

Subtraction of the expansions yields the velocity 

n{t)^^^m + St)^n{t-St)]+0{{St)^) . (15) 
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This Verlet algorithm [1,7,19,37,38] is manifestly time reversible: exchange of 
fi{t + St) and fi{t — St) yields the propagator for the time evolution going 
backward in time, and changes the sign of the velocity. 

A peculiarity of the Verlet algorithm is the fact that the updating of the 
velocities is one step behind: the velocities at time t can only be calculated, see 
Eq. H15|) . after the positions at time t^8t have become available. The updating of 
positions and velocities is synchronized in the "Velocity Verlet Algorithm" |40| . 
Using 

n{t + Si) = n{t) + 5tv, + ^('^O'/^W (16) 
together with the corresponding time-reversed equation 

n{t) = nit + St) - Stv,{t + St) + -L(^st)\nit + St) (17) 

Znii 

one obtains by adding these equations 

v,{t + St) = v,{t) + ^[Mt) + kt + St)] . (18) 

The propagation of the velocities hence requires that both the forces of the 
present and the future configurations are known. The second set of forces can 
be determined as soon as the coordinates at time t + St have been obtained. 
Note that both the Verlet and the Velocity Verlet Algorithms are completely 
equivalent: they produce identical trajectories. 

At this point it is of interest to ask what is the scale for the time step. For 
simplicity we consider liquid argon, the "fruit fly" for MD. The argon atoms are 
assumed to interact with a Lennard- Jones potential as considered in Eq. 
with parameters a ~ 3.4 A, s/ks ~ 120 K, and m « 6.6 x lO^^'^g. Rescaling 
coordinates by writing r * = r/a, we obtain 



f*{t + St)^2r*{t)-r*{t-St)-{St)' __ ^ ^ [(r*.)"" " (^I.rVS]. 

l'^^.?-! 

(19) 

From this equation we see that the natural time unit tq is 

/ ma'' \ 
'^-[-^) ' 

which is roughly tq « 3.1 x 10~^^ s for Argon. In order to keep numerical 
integration errors small, we need to choose the time step St so small that the 
third term on the rhs of Eq. p9|l is significantly smaller that the first and 
second term. For the Lennard- Jones system this usually means St* = 0.03, 
which translates into real time units St « 10~^^ s for the time step. So even 
with a million time steps (rescaled time t* = t/rQ = 30000) one only reaches a 
real time of about 10 ns. 



10 



It is also useful to go beyond this "pedestrian approach" to classical me- 
chanics, discussing more formally the time evolution in phase space, combin- 
ing the Cartesian coordinates of all the particles {X) and their momenta P = 
{pi,P2,--- ,Pn) into a point in a 2A^(i-dimensional space, T = {X,P). Liou- 
ville's theorem says that the flow in phase space has the property that it is 
incompressible, phase space volume is conserved. One can formally write for 
any observable A 

A(r,i) = c/(i)A(r,o), (21) 

where the propagator U{t) is a unitary operator, 

C/(i) =exp(i£t), (22) 

C being the Liouville operator. Now the Verlet algorithm has the desirable 
feature that (to order {StY) unitarity is preserved. This fact may to some extent 
explain why the Verlet algorithm is such a particularly stable algorithm over 
very long times (for a detailed discussion of the stability of the Verlet algorithm 
see Ref. ^). 

2.2 Estimation of Intensive Thermodynamic Variables; the 

NVT Ensemble 

If the energy E were rigorously conserved by the algorithm it would realize states 
distributed according to the microcanonical (NVE) ensemble of statistical me- 
chanics. Of course due to integration errors the energy is not strictly conserved 
and hence one has to make sure that the chosen time step is small enough, 
as will be discussed below - let us disregard this problem for the moment and 
assume that the microcanonical ensemble is realized perfectly. How do we then 
estimate the thermodynamic variables of interest, such as temperature T and 
pressure pi 

One can make use of the fact that in classical statistical mechanics the 
distributions of positions and velocities factorize; one simply has a Maxwell- 
Boltzmann distribution for the velocities and thus T can be inferred from the 
kinetic energy of the particles. Defining a "kinetic temperature" T as follows 
(for d = 3) 

2 1 ^ 

i—l 

the desired estimate of the temperature of the system is then computed as the 
average T = (T) . It must be emphasized that for finite N in the microcanonical 
ensemble there occur fluctuations in temperature, described by |42| 

(T)2 -3n[' 2cJ' ^ ^ 
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where is the specific heat of the system (at constant volume). Considering 
that one always works with a finite length of the MD run, these fluctuations 
cause a statistical error in the estimation of the temperature. 

The pressure p can be estimated using the virial Theorem. The dynamical 
variable whose average yields the pressure is defined as 




3N \ 



(25) 



such that 



N 

V ' 3V 

1=1 

Now we return to the problem that the energy is not strictly conserved by 
the algorithm described so far. The original recipe used in early MD work was 
to rescale the velocities from time to time such that one corrects for the drift in 
the energy. Alternatively, one can draw randomly velocities from time to time 
from a Maxwell Boltzmann distribution. Both methods have the disadvantage 
that dynamic correlations arc somewhat disturbed, and one realizes in this way 
neither a microcanonical nor a canonical ensemble. In the canonical (NVT) 
ensemble, temperature is a given variable and therefore not fluctuating at all, 
and instead of the fluctuations of T {Eq. H24|) l- one encounters fluctuations of 
the total energy E, described by another fluctuation relation: 

NCJkB = il/kBT)^{n^),,., - {H)l^^^) . (27) 

Here H is the Hamiltonian of the system. With methods in which one rescales 
the velocities, both temperature and energy are fluctuating, and neither Eq. H24|l 
nor Eq. ^ hold. 

However, one can introduce a MD method that does reproduce the canonical 
(NVT) ensemble exactly: one has to extend the Lagrangian of the system by 
a variable, representing the thermostat which has a flctitious mass Q. This 
yields the Nose-Hoover algorithm (43n44| . Newton's equations of motion are 
then extended by a "friction" term, 

fi/rrii- (fi . (28) 

The friction coefficient C(t) fluctuates in time around zero, according to the 
equation 



N 



(29) 



Thus, Q{t) responds to the imbalance between the instantaneous kinetic energy 
and the intended canonical average (remember ^ mi{vf) = SNksT). Although 
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the total energy is no longer conserved, one can identify a conserved energy-like 
quantity W, 



n' =H + -QC^ + SNksT / C{t')dt' 



(30) 



for which one can show that dTi' /dt ~ . 

Also in this case it is true that dynamical correlation functions between 
observables such as {A{Q)A{t)) are not precisely identical to the microcanon- 
ical ones. However, in many cases of practical interest (e.g., dense systems 
of flexible bead-spring-type chain molecules representing a glass-forming poly- 
mer melt IIS]) the difference between the result for {A{0)A{t)) from a strictly 
microcanonical run and a run using this Nose-Hoover thermostat is negligibly 
small pS] . 

It is also possible to realize the isothermal-isobaric (NpT) ensemble, where 
the pressure is given and rather the volume fluctuates, by coupling to the so- 
called "Andersen barostat" 0^1 • We shall not describe this here, but refer the 
reader to the literature [TirTllS^ for details. 

2.3 Diffusion, Hydrodynamic Slowing Down, and the Ein- 
stein Relation 

As a final point of this discussion of technical aspects, let us discuss how one 
extracts diffusion constants and other transport coefficients from MD simu- 
lations. We consider first the phenomenological description of diffusion in a 
single-component system in terms of Fick's law since this will allow us to dis- 
cuss concepts such as hydrodynamic slowing down, the Einstein relation for 
the diffusion coefficient D in terms of mean square displacements of particles, 
and the Green-Kubo formula providing a link with the velocity autocorrelation 

function umniiiTi. 

Fick's law states that there is a current of particles caused if there is a 
(particle) density gradient (or a gradient in the concentration) and this current 
acts to reduce the gradient. For small and slow density variations this is a linear 
relation. 



where D is the diffusion constant. If we combine this (phenomenological!) con- 
stitutive equation of irreversible thermodynamics with the (exact!) continuity 
equation which expresses the fact that the total particle number N is conserved, 



(31) 



dpir,t)/dt + \/-j{f,t)^0 
we get the well-known diffusion equation. 



(32) 



dp{r,t)/dt = D\/^p{r,t) 



(33) 
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This equation is easily solved both in real space and in reciprocal space. In 
reciprocal space a simple exponential relaxation of the Fourier components of 
the density fluctuations Spj^{t) results, 

6p{f, t) = p(f, t) - (p) = j Sp{k, t) exp(i k ■ r)dk , (34) 

since 

^ 5p(^k, t) = -Dk^Spik, t) , Sp{k, t) = 5p{k, 0) eM^Dkh). (35) 
We recognize that the associate relaxation time diverges as /c — * 0, 

= {Dk^)-\ (36) 

Therefore we see that the dynamic correlation function of density fluctuations 
for long wavelengths in a diffusive system decays very slowly, 

S{k, t) = {6p{~k, 0)Sp{k, t)) = S{k) exp{-Dkh) , (37) 

where S{k) is the static structure factor from Eq. S{k) = {6p{—k, 0)Sp{k, 0)). 
Eqs. 136|l . (|37|l demonstrate the so-called "hydrodynamic slowing down" |48| . 
On the one hand, analysis of the intermediate scattering function as given by 
Eq. H37|l allows to extract the diffusion constant. On the other hand, this hy- 
drodynamic slowing down is a difhculty for equilibration in the NVT ensemble, 
both for MD and for MC simulations. Remember that for the simulations of 
fluids, it is often convenient to use as an initial state a regular arrangement 
of the particles in the box, by simply putting them on the sites of the crystal 
lattice. This implies that S{k ^ 0) = in the initial state. If the simulated 
volume is large, so that the smallest wavevector fcmin = Stt/L is small, it will 
take a very long time until iS'(A;min) reaches its equilibrium value. In judging 
whether or not full equilibrium has been achieved, one hence cannot rely on 
the recipe advocated in some of the early simulation literature to check whether 
the internal energy has reached its equilibrium value. This was okay for the 
early work, where only 64 or 256 particles were simulated, and hence k^{^ was 
not small. Being interested in the simulation of much larger systems today, 
the issue of equilibrating long wavelength fluctuations properly needs careful 
consideration |42}. 

It is also useful to examine the diffusion equation in real space since it readily 
allows to derive the Einstein relation for the mean square displacement of a 
diffusing particle. Let us take as an initial condition for the solution of Eq. H33|l 
a J— function, p(f, t = 0) = 5(r — rg) {note that the normalization J p{f, t)dr = 1 
means that p{r^ t) can be interpreted physically as the conditional probability 
to find a particle at f after a time t provided this particle was ai f — fo at time 
t = 0}. Now the solution of Eq. (|33|l is simply a Gaussian distribution, 

p(r,i) = (47ri?t)-'^/2 exp[-(f-ro)V(4i?0] , (38) 
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d denoting once more the dimension of space. The squared half width of this 
distribution increases linear with time i, and so does the mean square displace- 
ment: 

{[5r{t)f) = {[r{t)-raf)^2dDt , t ^ oo. (39) 

We have added here the restriction to consider large times {t oo) since then 
this Einstein relation, Eq. 1)39(1 . holds quite generally, while typically the simple 
diffusion equation, Eq. H33|l . will not hold on small length scales (of the order 
of a few atomic diameters or less) and the corresponding short times. This fact 
will be evident with the examples discussed later. In fact, Eq. ijS^ is routinely 
used in simulations to compute the (self-)diffusion constants. 

2.4 Green Kubo Relations and Transport Coefficients 

It is also interesting to note that there is a relation between the self-diffusion 
constant and the velocity autocorrelation function of a diffusing particle (1^ . 
Writing 

t 

r{t)-ro= j vlt')dt' , (40) 



the mean square displacement can be expressed as follows: 



{[r{t) - fof) ^ dt' dt"{v{t") ■ v{t')) ^d dt' dt"Z{t" - 1'), (41) 



with Z{t" — t') = {va{t")va{t')). In this notation, we imply that in equilibrium 
translation invariance holds, so Z depends only on the time difference t" — t' , 
not on the two times t' ^t" separately. Va is one of the d Cartesian components 
of V. Rearranging the domain of integration in Eq. H41|l readily yields 

t 

{[r{t) - fo]2) = 2dt / (1 - s/t)Z{s)ds — * 2dDt , (42) 

J t^oo 


with 



J Z{s)ds = j {va{Q)vc{t))dt. (43) 



This result that the self-diffusion constant is the time integral of the velocity 
autocorrelation function is a special case of a "Green-Kubo relation" An- 
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other example of such a relation is that the shear viscosity rj is related to time 
correlations of the off-diagonal components of the pressure tensor axy. 



1 



' VksT 

where 



dt{<Jxy{0)axy{t)) , (44) 



^ f 1 1 

fy being the y-component of the force with which particles i,j interact, fy = 
■-dU{rij)/dy, cf. Eq. rij = fi — r}. Similarly, for the thermal conductivity 
At we need to consider correlations of the energy current density [T^I17| 

oo 

- ! dt{fMm) , (46) 







dt ' 

Finally, the electrical conductivity (Tei can be obtained from correlations of 
the current density Jf of the electrical charges qi |19II47| 

oc 

^ci = J dt{J^\0)J^\t)) , (48) 



N 
i=l 

All these relations are in fact useful for MD simulations. 



3 Application to the Example of molten Silica 

As an example for the type of results that we can get using the Einstein 
and Green-Kubo relations that were just discussed, we show in Fig. |21 an Ar- 
rhenius plot for the viscosity of molten Si02 ISOj, and compare the MD re- 
sults with experimental data (FT]. "Arrhenius plot" means, the logarithm of 
the viscosity is plotted vs. inverse temperature, since then an Arrhenius law 
[ri(T) cx exp[i?A/(fcs7')], where Ea is an "activation energy"] shows up as a 
straight line in the plot, with a slope E^- The experimental data included 
here |51| is indeed compatible with such a law. However, the extremely large 
values of the viscosity mean that the relaxation time of this melt is in the 
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microsecond or even millisecond range. Therefore the MD results cannot be ob- 
tained in the same temperature range where the experimental data were taken, 
but only at considerably higher temperatures, where experiments are no longer 
possible. Nevertheless, these simulation results (together with results for many 
other correlation functions [SUl that will not be shown here) do have a great 
theoretical relevance, since they show that even Si02 has a regime of tempera- 
tures where the so-called mode coupling theory of the glass transition [52] can 
be applied. Since this regime of temperatures is inaccessible with laboratory 
experiments, the hypothesis was advanced that one needs several classes for 
glass-forming fluids with a fundamentally different dynamical behavior. Fig. |31 
shows that "computer experiments" can complement laboratory experiments by 
providing results for physical quantities in a parameter range which can not be 
studied in the laboratory. Another interesting aspect of the data shown in Fig.O 
is that the Stokes-Einstein relation, linking the temperature dependence of the 
viscosity to that of the self-diffusion constant(s), does not hold. 

Simulating real glass-forming fluids such as Si02 is difficult for various rea- 
sons: a great problem is again the choice of a suitable force field. In Fig. O the 
so-called BKS-potential jSH was used, which has been proposed on the basis of 
quantum-chemical calculations. It contains Coulomb-like interactions, but with 
effective charges qi{i S Si, O } rather than the true ionic charges, and a short 
range Buckingham potential, 

2 

Uin,) = + A, eM-Brjn,) - C,,lr% , (50) 

where e is the elementary charge, go — ~l-2, q<s,\ — +2.4, and the constants yl^, 
Bij and can be found in the original reference jSS^ . It is somewhat surprising 
that a clever choice of these phenomenological constants allows to describe the 
directional covalcnt bonding (typically a Si atom is in the center of a tetrahe- 
dron, with the oxygens at the corners), although one uses just pair potentials. 
Nevertheless, the simulations are still technically very difficult: the long range 
Coulomb interaction necessitates the use of the time-consuming Ewald summa- 
tion techniques ^121131 \ the scale for the potential is in the eV energy range and 
varies rather rapidly with distance. Therefore one needs to use a rather small 
MD time step, namely 8t ~ 1.6 fs. 

If one produces amorphous glassy structures by a MD simulation with such 
a model, the most plausible procedure is the same as the one used in the glass 
factory. One starts at a very high temperature Tq, and then one cools the system 
gradually down at constant pressure, so the temperature Tify varies with time, 
e.g. T(t) = Tq — 7t, 7 being the cooling rate. While the structures obtained in 
this way look qualitatively reasonable [SJ , one must be aware of the fact that 
the cooling rates applied in the simulation are extremely large (7 = 10^^ K/s or 
even higher), which means that they are at least a factor 10^^ larger than those 
used in the experiment. Thus, one must expect that the quantitative details 
of the results (including also the density of the material) will depend on the 
cooling rate distinctly, and this is what actually is found (54|. 
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Better results are obtained if one proceeds in a slightly different way, where 
one no longer tries to predict the density of the model system from the simula- 
tion, but rather fixes it to the experimental value. The system then is carefully 
equilibrated at a "moderately high" temperature before one cools it down to 
the desired temperature. "Moderately high" for pure Si02 means T = 2750 
K - then equilibration requires a run which is already 20 ns long. To avoid 
finite size effects, a box of around 48 A linear dimension needs to be taken, 
containing about 8000 atoms. This calculation requires substantial CPU re- 
sources (it was done [501 at a CRAY-T3E parallel supercomputer, performing 
force parallelization) so it is presently not easy to do much better. 

In any case, with this procedure one does obtain the static structure factor of 
silica glass at room temperature in very good agreement with experiment. This 
can be inferred, e.g, from Fig. 0] where we show the static structure factor as 
measured in a neutron-scattering experiment and which is the weighted sum of 
the three partial static structure factors [1915015 1| . Since there are no adjustable 
parameters whatsoever involved in this comparison, this agreement is significant. 

Fig. |S1 shows now an Arrhenius plot for the self-diffusion constants, using 
scales such that the experimental data |5S1|S7] can be included together with 
the simulation results. From this graph it becomes obvious that in principle 
one would like that the simulation spans over 16 decades in dynamics, a task 
clearly impossible for MD. So there is again a gap between the range where 
simulations in thermal equilibrium are possible and the temperature range where 
the experimental data can be taken. But it is gratifying that the simulation can 
predict the activation energies for the self-diffusion constants almost correctly. 

Now the real strength of the MD simulations is that one can record in a single 
calculation a great variety of different static and dynamic properties simulta- 
neously. For instances, one can study the time dependence of autocorrelation 
functions of the generalized temperature fluctuation 5T^{t), which is the Fourier 
transform of the local "temperature" T{r,t), 



One sees (Fig. that for large wavelength q this correlator is rapidly 

relaxing but again one notes a slowing down at small q. This is another exam- 
ple for "hydrodynamic slowing down". Of course, the considerable statistical 
scatter of the results shown in Fig. is somewhat disturbing. But it has to be 
mentioned that these results are already based on averages over 100-200 runs, 
corresponding to an effort of several years if one would run the problem on a 
single CPU, and hence it is currently difficult to do better. 
Defining now from the time-displaced correlation function. 




(51) 



(52) 



a relaxation time r by the condition 



^TT{q,t = t) = ^TT{q,t = 0)/10 



(53) 
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one obtains the data shown in Fig. EId. Here a plot of rq^ vs. is presented, 
since theory predicts jSH] 



r = 




(54) 



where p = N/V is the density, Cp the heat capacity, and At is the thermal 
conductivity. Since p is given and Cp can be estimated independently, one ob- 
tains a rough estimate of the thermal conductivity At from these data (Fig. EJd) . 
This estimate is in reasonable qualitative agreement with experiments, which 
indicate that for T > 1000 K At is only weakly dependent on temperature and 
in the range 2 < At < 3 W/Km [55). 

As a last example, the temperature dependence of the longitudinal (cf) and 
the transverse (ct) sound velocity is shown in Fig. [7| 59j. These results were 
obtained from an analysis of the corresponding time-displaced current-current 
correlation functions which were Fourier transformed into frequency space. Un- 
damped propagation of sound manifests itself by 5— functions 6{lu — ujq) with 
ojq = ci^tq for g — > in these correlators (SHI- Of course, in a liquid for g ^ 
no static shear can be maintained and hence in that case only q is well-defined. 
However, for not too small q both longitudinal and transverse correlators show 
broad peaks, and the positions of these peaks are shown in Fig. Again one 
notes very nice agreement with corresponding experimental data |()0j . 



4 A Brief Introduction to Non-Equilibrium Molec- 
ular Dynamics (NEMD) 



Already in Eq. H31|l . we have considered the situation that the density in the 
system may deviate from its constant equilibrium value, and we have postulated 
a linear relation between the particle density current and the gradient of the local 
density Vp(f, t). Similarly, we also can assume that there is no longer complete 
thermal equilibrium, as far as the temperature is concerned, but only "local 
equilibrium" : we may still assume a Maxwell-Boltzmann velocity distribution 
of the particles but with a temperature T{r, t) that slowly varies in space. Thus, 
a current of energy, i.e. heat, is created, the coefficient between the temperature 
gradient and the energy current density being the thermal conductivity At, 
yielding Fourier's law of heat conduction. 



Now in reality the situation is not so simple, since energy density and particle 
number density are coupled variables. Each gradient therefore produces also a 
current of the other variable. So we have to generalize the set of flow equations 
Eq. (|31|l . (|55|l to a matrix form, involving the Onsager coefficients A^^, 



jQ = -AT[vr(f,i)]. 



(55) 




(56) 
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with being the chemical potential of particle species i (the generalization of 
Eqs. H56() . H57() to systems containing several different species of particles then 
is obvious). Also, we have written the set of constitutive equations entirely 
in terms of gradients of intensive variables, rather than using the density of 
an extensive variable as in Eq. H31|) . Of course, using thermodynamics, the 
corresponding relation is easily established, since with Eq. H31|) on obtains: 

"^^^ ^^^'^ ' Vp. ^D,^p^(^±l\ . (58) 



ksT \dpi J ksT ksT \dpi 

In NEMD simulations one usually sets up a stationary gradient by suitable 
boundary conditions creating a stationary current through the system. In the 
case of a flow of particles, this clearly is compatible with periodic boundary 
conditions: particles leaving the simulation box at the right simultaneously 
reenter at the left. For instance, we may confine a fluid between two parallel 
plates and let a constant force act on the particles pnil21ll(jllfi2lf):^j but schemes 
where one avoids external walls altogether also are possible [lOlEllS]. In the 
latter case, one can obtain a homogeneous flow, while in the former case the 
velocity profile is inhomogeneous, and also the structure of the fluid close to the 
external walls is modified. 

Now an important aspect of all such stationary flows is that entropy is pro- 
duced: the entropy production per unit time is 



ds 



V 



knT 



+Ai. 



ksT 

-\ 2 



keT 



(59) 



Although this entropy production per unit time is small for small gradients, 
nevertheless the heat that is produced makes a strictly microcanonical simula- 
tion impossible: the system would heat up steadily. Thus one always has to use 
a thermostat. 

As an example for such applications of NEMD, we now briefly review some 
results jn21 obtained for a bead-spring model of short polymer chains confined 
between two parallel plates. The interaction between the beads is a simple 
Lennard- Jones (LJ) potential of the form as written in Eq. ® {truncated at 
Tc = 2.24(7 and shifted so that U{r = Vc) = 0, as usual}. The "spring poten- 
tial" between neighboring effective monomeric units is taken to be the finitely 
extensible nonlinear elastic (FENE) potential. 



C/fene(^) - -15i?2ln[l - [l/Rof] 



(60) 
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choosing "Lennard- Jones units" for length, energy and temperature as usual, 
(T = l,£ = l, ^5 = 1 {and also the mass of the effective monomers m = 1, 
which fixes the unit of time as well, cf. Eq. H2(J|) T. Choosing Rq — 1.5 one 
finds that the total bond potential (LJ + FENE) has its minimum at Iq « 0.96, 
while the LJ potential has its minimum at tq = 2^/^ « 1.13. This competition 
between these two length scales prevents crystallization of this melt, and hence 
is the physical reason that this model is a very good model for a glass-forming 
polymer melt |lHlinS]. Note that unlike the atomistic model (34j|35j discussed 
in Sec. 1.3 of the present article, neither a bond angle potential nor a torsional 
potential is included in this case. In fact, one envisages that the "effective 
monomers" are formed by combining a few subsequent chemical monomers along 
the backbone of the chain molecule into one such effective unit of a "coarse- 
grained" model HllllEH]. 

For a chain length iV = 10 a box of linear dimensions Lx x Ly x D, with 
Lx — Ly — 10.05, D — 20, containing 200 chains in the system is a good 
choice |62| . We used periodic boundary conditions in x and y directions, while 
in the remaining ^-direction one places walls formed by atoms on a triangular 
lattice, at equilibrium positions r^^oq with Zi^cq — iD/2. These wall atoms are 
bound to their equilibrium positions with a harmonic potential, C/harm(^i) = 
^Khin - r^,oq)^ with Kh = 100. 

In order to represent the physical effect of the atoms in the (massive) walls 
other than those in the first layers of the wall, one coarse-grains them into an 
effectively repulsive background potential U^/ediiz) — e{a/z)^, where z = \zaion — 
-Zwaiil, Zaion being the z-coordinate of the considered monomer, and Zwaii = 
zt{a+D/2) the effective positions of the second layer on top of the walls. In order 
to create a flow in this geometry, a force F'^ acts in the -|-x-direction on each 
monomer, and one also needs to specify the interaction between the monomers 
and the atoms in the top layers of the walls, which determine the "hydrodynamic 
boundary conditions" of the resulting Poisseuille flow pUll^l^l^ . 

In Ref. 1^2] this interaction was chosen of the same LJ form as the interaction 
between non-bonded effective monomers, but with different parameters, CTwm = 
2~^/^, £wm = 2. In this way one creates a "stick" boundary condition for the 
flow, while for awm = 1, £win — 1 one would have a strong partial slip, and the 
estimation of transport coefficients from the velocity and temperature profiles 
would be more difficult. 

Fig.lHlshows typical results for the density and velocity profile across the slit. 
One can recognize a very pronounced "layering" phenomenon near both walls, 
i.e. p{z) exhibits very strong oscillations, but these oscillations are only weakly 
dependent on temperature, and furthermore in the center region of the thin 
polymer film (—5 < z < -f 5) the density profile is essentially flat, exhibiting the 
behavior of the bulk. Using the radial distribution function g{r), with a vector 
r parallel to the walls, it has been checked that the system has still a fluid like 
structure (and not crystallized, for instance) near the walls where the strong 
layering phenomena occurs. 

The central region (—5 < z < +5) then is used to fit the Poisseuille-type 
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flow profile (well-known from hydrodynamics) to the data, i.e. 



Ux{z) 



'wall 



2Zwall<5)/»7 



(61) 



where po = p(z = 0) is the bulk density, and 6 has the physical interpretation 
of a slip length {note dux{z)/dz \ ~ Ux{z = z^a\\)/5}. Fig. 13) shows that 



excellent fits to Eq. Ht)l|) are in fact obtained, thus yielding an estimate of the 
viscosity rj. 

At this point a comment on how the thermostating is done is in order. The 
most natural way (used in the laboratory experiments) would be to enforce 
isothermal conditions through the walls, keeping them at constant temperature 
by coupling the wall atoms to a thermostat. The heat created in the center of 
the polymer film through the viscous flow then leads to a stationary temperature 
profile, and along the resulting temperature gradients the heat is transported 
towards the walls. Of course, if the viscous flow is too fast (when the force 
amplitude F"^ is too large) one enters a regime of very nonlinear response, with 
strong temperature variations and large gradients, and then descriptions based 
on linear irreversible thermodynamics {Eqs. H51|) ~ (|61|l } are not really ap- 

plicable in a strict sense. As an example, Fig. |5| shows temperature profiles 
across the film observed in a simulation of this type. One sees that for F"^ > 0.1 
the temperature in the center of the film is strongly enhanced, and one needs to 
choose < 0.05 to stay within the linear response regime. Of course, the study 
of nonlinear phenomena far from equilibrium (such as "shear thinning" j21| , the 
decrease of the effective viscosity r^eS with increasing shear rate) may be of 
interest in itself, but this is out of consideration here. 

Fig. 121 also demonstrates that the temperature profiles can also be fitted 
to the theoretical profiles resulting from the solutions of the phenomenological 
hydrodynamic equations j65ll67j 



In the temperature range of Fig. one obtains 3 < At < 4.4 as a result for 
the thermal conductivity. We also note that the data shown in Fig. [HI do not 
give any evidence for an additional term not included in Fourier's law, Eq. (|55|l . 
namely a strain rate coupling that would lead to an additional quadratic term 
{oc -{2z/Df} in Eq. ^ 61 . 

Since in a glass forming fluid the viscosity is strongly temperature dependent, 
already a small variation of T{z) across the film creates problems for the use of 
Eq. (|62|) at low temperature. It thus is preferable to apply the thermostating 
algorithm not only to the wall atoms, but also to the fluid atoms, noting that in 
the presence of the flow the temperature is defined as T = m{{v — (t/))^)/(3fcs) 
instead of T = (T) (Eg. f 1^ 1. Since here (w) = {ux{z),Q,Q), one has during 
the equilibration to determine the profile Ux{z) self-consistently |62| . 

Fig. 1 101 shows the result for the shear viscosity 77 and the lateral self-diffusion 
constant D±^ (obtained from mean square displacements in the y-direction per- 




(62) 
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pendicular to the flow). One sees that one can obtain rather precise results for 
r] and monitor the increase of rj over about two decades. Both 77 and D± can 
be fitted well by the Vogel-Fulcher-Tammann (VFT) relation, 



7KT) = r7(w) exp[£;^/(T ~ To)] , = D{^) eM-Eo/iT - Tq)] , (63) 

where the VFT temperature Tq is the same for both quantities, while the effec- 
tive activation energies i?,,, Ed are slightly different. This difference is consistent 
with a direct examination of the ratio I = kgT/ {ATrD±ri), which should be a 
characteristic constant length if the Stokes-Einstein relation holds (cf. Fig.O: 
similarly as for Si02, one finds a systematic decrease of this ratio as the tem- 
perature gets lower [5^ . 



5 Concluding Comments 

In this brief review, we have attempted to convey to the reader the flavor of MD 
simulations, addressing both the estimation of equilibrium properties such as 
static structure factors (Fig. 0)) or density profiles near confining walls (Fig. |Sli) 
and of dynamic properties, such as single chain intermediate dynamic structure 
factors in polymer melts (Figs. ^ 13 1 diffusion constants in Si02 (Fig. [3) and 
glass forming polymer melts (Fig. llOb ). and collective transport coefficients 
such as viscosity (Fier. [TOkl. thermal conductivity, etc. Both type of MD 
simulations, dealing with systems in thermal equilibrium and NEMD, has been 
discussed (Sec. 4). 

Already with respect to static properties, one has the choice between differ- 
ent ensembles (microcanonical NVE ensemble, or NVT and NpT ensembles, 
realized with suitable thermostats and barostats). It depends on the questions 
which one likes to answer which ensemble is more appropriate. Similarly, we 
have demonstrated that transport coefficients (such as viscosity, thermal con- 
ductivity, etc.) can be estimated from thermal equilibrium simulations (via the 
Green-Kubo jBS] relations of Sec. 2.4) or via NEMD work. There are many 
variants how one can proceed |2niEIlEllEllEniEIlE2E3IMIEZI , and clearly 
there is no full consensus in the literature about the "pros" and "cons" of the 
various approaches yet. Instead the subject is still a matter of current research. 
As an example, Fig. II II presents an example, again Si02, where viscosities were 
estimated from a NEMD approach in complete analogy with Fig. and the 
data are compared to the Green-Kubo approach used in Fig.|2| At least within 
the error bars noted in Fig. |31 there is very good mutual agreement, and this 
is reassuring because it shows that the various systematic errors of the simula- 
tions and their analysis are under control. We have emphasized that one must 
be aware of systematic errors due to the finite size of the simulation box, due 
to incomplete equilibration, in particular in the case of long wavelength fluctu- 
ations, due to inaccuracies of the algorithm related to the size of the time step, 
etc. 
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Finally, we emphasize that we have used illustrative examples which were 
taken from the research of the authors for the sake of simplicity only — similarly 
valuable research on related problems is available in the literature from many 
other groups, as is well documented (see e.g. Refs. P1I^3I7I19II2UII21) '). We hope 
that the present article "wets the appetite" of the reader to study this extensive 
literature. 
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Figure 1: Intermediate dynamic structure factor versus scaled time as obtained 
from a neutron spin echo experiment (symbols) and computer simulation results 
(full curves). The time axis is scaled by the respective center of mass diffusion 
coefficients Dn. The ordinate axis is scaled by the static (equal-time) single 
chain structure factor S{q,Q). The different symbols represent the different 
values of q, in units of A^^, as explained in the legend. From Paul |35|. 
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Figure 2: Normalized coherent intermediate scattering function for C100H202 
plotted versus time, for the wavelength q (in A~^) as indicated in the legend. 
The full curves are the Rouse model prediction, the symbols denote the simu- 
lation results. From Paul \6b\ . 
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Figure 3: Molecular Dynamics results (filled squares) for the viscosity of the 
BKS model for molten Si02 plotted vs. inverse temperature. The dashed 
straight line indicates an Arrhenius fit with an activation energy Ea — 5.19 eV. 
The experimental data |51| (open circles) are compatible with this value but 
would suggest a slightly larger preexponential factor. Note that for Si02 the 
analysis of other correlation functions at very high temperature suggests a crit- 
ical temperature = 3330 K of mode coupling theory, therefore r]{T > Tc) 
deviates from the Arrhenius law. The inset shows that the Stokes- Einstein rela- 
tion, kBT/{rjD) — const, where D are the Si or O self-diffusion constants, does 
not hold in the regime of temperatures studied. From Horbach and Kob |50|. 
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Figure 4: Static neutron structure factor of Si02 at room temperature (T = 
300 K), plotted versus wavelength q. The full curve is the MD simulation |5U|. 
using the experimental scattering lengths for Si and O atoms, while the symbols 
are the neutron scattering data of Price and Carpenter From Horbach and 

Kob isni. 
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Figure 5: Plot of the self-difFusion constant D of silicon atoms (Si) and oxygen 
atoms (O) in molten Si02 as a function of inverse temperature. The symbols in 
the upper left part are the results from MD simulations and the data in the lower 
right part stems from experiments j56[l57) . The thin straight lines show simple 
Arrhcnius behavior {D oc exp[— i?A/(fcs2^)]) with an activation energy E^, as 
indicated in the figure. The vertical broken lines indicate the experimental 
glass transition temperatures, Tg = 1450 K, as well as values for Tg that one 
obtains if one extrapolates the data from the simulations to low temperatures 
and then estimates Tg from the experimental value of the O diffusion constant 
{DoiT = T|"") = 10-16 cmVs =^ T|'™ ^ 1381 K ) or the Si diffusion constant, 
respectively (DsiiT ^ = S.IO^^^ cmVs =^ T|™ = 1303 K). From 

Horbach and Kob 



32 



0.5 

C\i 

50.4 

I- 

^ 0.3 
0.2 
0.1 
0.0 

0.04 









T=3000K 


\ Vi 

\ 

\ \\\ 
\ \\ 
\ \\ 
\ """^ 


T=3580K 

k T=4700K 


\ \ 

\ V 

'v\q=0.37A"' 

\V 


\ 

\ r 


\q=0.84A:\* 


q=1 .75A"' 


y V 



0.01 



0.00 



10" 



t[ps] 




0.0 



1.0 



2.0 



q [A] 



Figure 6: a) Time-dependence of the autocorrelation function ^TT{<l,t) /T^ of 
the generahzed temperature fluctuation STq{t) for Si02, showing data for three 
temperatures T and three choices of wavevector q as indicated. AU data refer 
to systems containing 224 oxygen and 112 sihcon atoms, b) Plot of wavevector 
squared times the relaxation time t, as defined in Eq. H53|l . versus . From 
the extrapolated intercept for q^ 0, the estimate At ~ 2.4 W/(Km) for the 
thermal conductivity At is extracted. From Scheidler et al. |58| . 
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Figure 7: Longitudinal (ce) and transverse (ct) sound velocities of Si02 plotted 
vs. temperature. These quantities were determined from the frequency posi- 
tions of the maxima of the corresponding longitudinal and transverse current 
correlation functions at q — 0.13 A^^ (open circles) and q = 0.18 (filled 
squares) . Also included are the experimental data of Polian et al. ^Uj which are 
multiplied with the factor (2.2/2.37)^/^ since the simulation was done at a den- 
sity of psim = 2.37 g/cm^ while the experiment was done for poxp = 2.2 g/cm^. 
From Horbach et al. |59j . 
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Figure 8: a) Monomer density profiles of a polymer melt confined in a slit. 
The bulk density is p = p{z = 0) = 0.99 and the amplitude of the force in 
the x-direction is F'^ = 0.05. Three temperatures are included as indicated. 
Note that the coordinate origin for the z-axis is chosen in the center of the 
polymer film, b) Velocity profile Ux{z) for two temperatures (as indicated) and 
otherwise the same conditions as in panel a) . The parabolic curves indicate the 
fitted Poisseuille flow profiles, Eq. Ht)l|) . From Varnik and Binder (62| . 



35 



] — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — f 




J I I I I I I I I I I I I I I I I I I I L 



-10 -5 5 10 

Z 

Figure 9: A comparison of the temperature profiles resulting from NEMD 
simulations (symbols) with the theoretical prediction Eq. (|62|l (lines) for various 
values of i^*^ at a wall temperature Twaii = 1 and a bulk density p{z = 0) = 0.795. 
Contrary to the results shown in Figs. |H1 and EH here the fluid particles were 
not coupled to a heat bath but obeyed pure Newtonian dynamics, while the 
walls were thermostated. Note that in computing the local temperature at z 
the streaming velocity u{z) is subtracted from the instantaneous velocities of all 
particles in the interval [z, z + dz]. From Varnik and Binder |()2| . 
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Figure 10: a) Shear viscosity of the model polymer melt as a function of the tem- 
perature distance T — Tq from the Vogel-Fulcher-Tammann (VFT)-temperature 
To = 0.190 ± 0.005. Symbols denote the data obtained for constant density 
Po — 0.99. The curve is a fit to the VFT equation {Eq. Ht)3|) T. b) Same as a) 
but for the selfdiffusion constant D±^ (measured in the y-direction transverse to 
the flow). From Varnik and Binder (62] . 
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Figure 11: Velocity profile of molten Si02 between atomistic walls, for an 
acceleration ax — 9.6 A/ps^, and three temperatures as indicated. Points are 
NEMD data using a simulation box y< Ly x D, with — Ly ~ 23.0 A, 
D — 31.5 A, using altogether N = 1152 atoms. Curves are fits to Eq. (|61|l . The 
resulting viscosities are quoted in the figure, together with the corresponding 
Green-Kubo estimates t/gk from Ref. [SU]. From Horbach and Binder |63| . 
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